## To create and send PBS files to run natural waterbody simulations ###

setwd("~/final_natural_water_files/PBS_Files")

script.file<-'~/evo-dispersal/art_wbdies/Pilbara/Pilbara_w_artificial.R'

for (ss in 1:100) {
  
    
    	fid<-ss
	    ##create the sh file
	    zz = file(paste("Pilb_Art", fid,'.sh',sep=''),'w')
  	  cat('##################################\n',file=zz)
	    cat('#!/bin/sh\n',file=zz)
	    cat('cd $PBS_O_WORKDIR\n',file=zz)
  	  cat("R CMD BATCH --no-save --no-restore '--args File.ID=", fid, "' ", sep="", file=zz)
	    cat(script.file, " ", paste("Pilb_Art", fid,'.Rout',sep=''), "\n", sep="",file=zz)
	    cat('##################################\n',file=zz)
  	  close(zz)
			
    	#submit the job
    	if (ss==100) {system(paste("qsub -m e Pilb_Art", fid,".sh",sep="")); break}
	    system(paste("qsub -m n Pilb_Art", fid,".sh",sep=""))
  
}

########################
# And then concatenate the output
setwd("~/final_natural_water_files/Pilbara")
load(file="Natural and artificial waterbodies table.RData")
concat.times<-c()
concat.pops<-c()
for (ii in 1:100){
  load(file=paste("times colonised natural and artificial", ii, ".RData", sep=""))
  concat.times<-c(concat.times, times)
  concat.pops<-rbind(concat.pops, points.colonised)
}

# proportion that colonised
sum(concat.times<100)/1000
summary(concat.times)
sd(concat.times)

#create a histogram of times
png(file="Colonisation times.png", width=16, height=16, res=300, units="cm")
hist(concat.times[concat.times<100], bty="l", xlab="Colonisation time (years)", 
     ylab="Density", freq=F, col="grey80", main="", ylim=c(0, max(density(concat.times[concat.times<100])$y)))
lines(density(concat.times[concat.times<100]))
dev.off()





